PRO MAKE_TEMPERATURE_FIGURE, suff
  
  spawn,'./irdc_dist_model/data/zzz_bring_calculon.sh'
  
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  IF n_elements(suff) EQ 0 THEN suff = ''
  
  IF ~ FILE_TEST('./irdc_dist_model/data/case3_irdc_grid'+suff+'_AT.sav') THEN $
     CONSOLIDATE_FFORE_GRID,suff
  
  restore,'./irdc_dist_model/data/case3_irdc_grid'+suff+'_AT.sav',/ver
  
  
  ;; Delete this file if n != 32
  fn = FILE_SEARCH('./irdc_dist_model/data/case3_irdc_grid'+suff+$
                   '_block*.sav',count=nf)
  IF nf NE 32 THEN $
     spawn,'rm ./irdc_dist_model/data/case3_irdc_grid'+suff+'_AT.sav'
  IF n_elements(rb3) EQ 0 THEN $
     restore,'./irdc_dist_model/bgps_rb3'+suff+'.sav',/ver
  n = n_elements(rb3)
  
  pcolor = 'BLK4'
  
  myps,'./irdc_dist_model/analysis_plots/temp_analysis'+suff+'.eps',$
       xsize=10,ysize=4.5,/cmyk
  
  multiplot,[2,1],xgap=0.035
  
  ;;==================================================
  ;; Frist Panel -- Temperature Histogram
  
  ind = where(rb3.batt_td NE 0, nind)
  print,nind
  
  print,m4_stat(rb3[ind].l)
     
  bin = 1.
  plothist,rb3[ind].batt_td,xarr,yarr,charsize=1.0,thick=3,bin=bin,$
           xtit='Dust Temperature from Battersby et al. (2011)',$
           ytit='N per '+string(bin,format="(F0.1)")+' K bin'
  yfit = mpfitpeak(xarr,yarr,A,nterms=3)
  print,A
  xp = findgen(201) / 200. * (max(xarr) - min(xarr) + 2*bin) + min(xarr)-bin
  oplot,xp,gauss_1(xp,A),thick=6,color=cgColor(pcolor)
  plothist,rb3[ind].batt_td,xarr,yarr,charsize=1.0,thick=3,bin=bin,$
           xtit='Dust Temperature from Battersby et al. (2011)',$
           ytit='N per '+string(bin,format="(F0.1)")+' K bin',/over
  axis,xaxis=0,xthick=3,xtickformat='blank_axis',/xst
  
  al_legend,/top,/left,thick=6,color=pcolor,box=0,linestyle=0,$
            ['$\mu$ = '+string(A[1],format="(F0.1)")+' K'],linsize=0.6
  al_legend,/top,/right,box=0,['N = '+string(nind,format="(I0)")]
  
  multiplot,/doxaxis,/doyaxis
  
  ;;==================================================
  ;; Second Panel -- A-Td isoprobability contours
  
  plot,[0],[0],/nodata,xr=[0,0.5],yr=[10,30],xtit='Foreground "veil" (A)',$
       ytit='Dust Temperature [K]'
  
  vline,/horiz,25,color=cgColor('Blue Violet')
  vline,/horiz,20,color=cgColor('Blue Violet')
  vline,/horiz,15,color=cgColor('Blue Violet')
  vline,0.13,color=cgColor('Blue Violet')
  vline,0.33,color=cgColor('Blue Violet')
  
  
  n_A = 51L
  A = dindgen(n_A)/100.d
  n_Td = 81L
  Td = dindgen(n_Td)*0.25d + 10.d

  levels = min(grid) + [2.3, 6.2, 11.8, 20, 30, 40, 50, 100]
  colors = [replicate('Opposite',3),replicate('RED6',4),'Opposite']
  print,levels
  
  contour,grid,A,Td,levels=levels,/overplot,c_colors=cgColor(colors)
 
  print,'Min chi^2_red = ',min(grid)/n_elements(rb3)
  
  myps,/done,/mp
 
  
END
